# Project: Serbia Survey
# Authors: William O'Brochta and Patrick Cunha Silva
# Goal: Analysis on voters' own attitudes (H1)

# Load Data
load(file = here('VotersOwnAttitudes.RData'))

##########################################################
################# Main Models for H1 #####################
##########################################################

#------------------ Unified DV ------------------------#
# Model for H1, no controls
m1a <- lm(Y ~ cyrillic_second, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))
ub1a <- qt(0.975, df = df.residual(m1a)) * se1a['cyrillic_second'] + coef(m1a)['cyrillic_second']
lb1a <- qt(0.025, df = df.residual(m1a)) * se1a['cyrillic_second'] + coef(m1a)['cyrillic_second']

#------------------ better_world ------------------------#
# Model for H1, no controls
m1b <- lm(better_world ~ cyrillic_second, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))
ub1b <- qt(0.975, df = df.residual(m1b)) * se1b['cyrillic_second'] + coef(m1b)['cyrillic_second']
lb1b <- qt(0.025, df = df.residual(m1b)) * se1b['cyrillic_second'] + coef(m1b)['cyrillic_second']

#------------------ most_places ------------------------#
# Model for H1, no controls
m1c <- lm(most_places ~ cyrillic_second, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))
ub1c <- qt(0.975, df = df.residual(m1c)) * se1c['cyrillic_second'] + coef(m1c)['cyrillic_second']
lb1c <- qt(0.025, df = df.residual(m1c)) * se1c['cyrillic_second'] + coef(m1c)['cyrillic_second']

#------------------ influence ------------------------#
# Model for H1, no controls
m1d <- lm(influence ~ cyrillic_second, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))
ub1d <- qt(0.975, df = df.residual(m1d)) * se1d['cyrillic_second'] + coef(m1d)['cyrillic_second']
lb1d <- qt(0.025, df = df.residual(m1d)) * se1d['cyrillic_second'] + coef(m1d)['cyrillic_second']

# Plot (Figure Average Treatment Effect of the use of Cyrillic Alphabet in the Survey on Voters’ Attitudes towards Nationalism)
par(mar = c(5, 12, 2, 2))
plot(y = c(0.75, 2, 3.25, 4.5), x = c(coef(m1b)['cyrillic_second'], 
                                      coef(m1c)['cyrillic_second'], coef(m1d)['cyrillic_second'], coef(m1a)['cyrillic_second']), 
     xlim = c(-0.4, 0.4), ylim = c(0, 5),
     pch = 20, cex = 2.25, 
     xlab = 'ATE of Cyrillic Script (Survey)', ylab = '', cex.lab = 1.75, 
     axes = FALSE)
lines(y = c(0.75, 0.75), x = c(ub1b, lb1b), lwd = 1.5)
lines(y = c(2, 2), x = c(ub1c, lb1c), lwd = 1.5)
lines(y = c(3.25, 3.25), x = c(ub1d, lb1d), lwd = 1.5)
lines(y = c(4.5, 4.5), x = c(ub1a, lb1a), lwd = 1.5)
axis(1, las = 1, cex.axis = 1.5, at = c(-0.4, -0.2, 0, 0.2, 0.4)
     , labels = c(-0.4, -0.2, 0, 0.2, 0.4))
axis(2, labels = c('Better World\nif Equal to Serbia', ' Serbia is Better\nThan Most Places', 'More Serbian\nInfluence is Better', 'Nationalism\nIndex'), 
     at = c(0.75, 2, 3.25, 4.5), tick = FALSE, cex.axis = 1.5, las = 2)
abline(v = 0, lty = 2)


###############################################################################
################# SI A.1 - H1: Descriptive Stats - Part 1 #####################
###############################################################################

# This only includes the dependent variables. The remaining variables come from file #3.

# List of DVs
dvs <- c("Y",
            "better_world",
            "most_places", 
            "influence")

# Covariates names
var.names <- c("Nationalism Index (Self)", 
               "More Serbian influence is Better (Self)",
               "Serbia is Better than Most Places (Self)",
               "World would be Better if it was Equal to Serbia (Self)")

# Generate Table
stargazer(mydf[, dvs], covariate.labels = var.names, summary = TRUE)


#######################################################################
############ SI A.4: Use of Cyrillic and Latin Alphabets ##############
######################################################################

## Use of Cyrillic and Latin at Home
round(prop.table(table(latin = mydf$latin_home, cyrillic = mydf$cyrillic_home)) * 100, 2)

## Use of Cyrillic and Latin at Work
round(prop.table(table(latin = mydf$latin_work, cyrillic = mydf$cyrillic_work)) * 100, 2)

#########################################################################
################# SI A.5 - H1: Randomization Checks #####################
#########################################################################

# Create Dummies
mydf <- fastDummies::dummy_cols(mydf, select_columns = c('gender', 'ageCat', 'other_languages'))

# List of covariates
covariates <- c("cyrillic_start",
                "gender_woman",
                "gender_man",
                "gender_no_response",
                "college",
                "ageCat_18-35",
                "ageCat_36-49",
                "ageCat_50-64",
                "ageCat_65+",
                "married",
                "unemployed",
                "urban",
                "other_languages_0",
                "other_languages_1",
                "other_languages_2",
                "other_languages_>2",
                "dif_lc_home",
                "dif_lc_work",
                'cyrillicImportance',
                "serbian",
                "sns",
                "polKnowQ1",
                "polKnowQ2",
                "intPolitics")

# Names of Variables
var.names <- c('Alphabet Ad (Cyrillic)', 
               "Gender = Woman",
               "Gender = Man", 
               "Gender = Prefer not to share", 
               "Education = College", 
               "Age = 18-35", 
               "Age = 36-49", 
               "Age = 50-64", 
               "Age = 65+", 
               "Marital status = Married", 
               "Unemployed", 
               "Urban", 
               "# of languages spoken = 0",
               "# of languages spoken = 1", 
               "# of languages spoken = 2", 
               "# of languages spoken > 2", 
               "Use of Cyrillic-Latin (at home)", 
               "Use of Cyrillic-Latin (at work)", 
               "Importance of Cyrillic",
               "Serbian", 
               "SNS",
               "Knows the # of MPs", 
               "Knows the presidential term length", 
               "Interest in politics")

# Save pvalues
pvalues <- list()
for(i in 1:length(covariates)){pvalues[[i]] <- t.test(mydf[, covariates[i]] ~ mydf$cyrillic_second)$p.value}
unlist(pvalues)

# Create a plot
par(mar = c(6, 15, 1, 1))
plot(x = pvalues, y = 1:length(covariates), 
     xlim = c(0, 1), 
     ylim = c(1,length(covariates)),
     pch = 20, 
     xlab = 'p-value',
     ylab = '', 
     axes = FALSE, 
     cex = 1.75, 
     cex.lab = 1.5)
axis(1, at = seq(0, 1, by = 0.05))
axis(2, at = 1:length(covariates), labels = var.names, las = 1)
abline(v = 0.01, lty = 3)
abline(v = 0.05, lty = 3)

#########################################################
############ SI A.6 - H1: Complete Results ##############
#########################################################

# Table for SI
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     style = 'apsr', type = 'text', omit = 'Constant')

#######################################################################
############ SI A.7 - H1: Descriptive Figures (Outcomes) ##############
#######################################################################


######### “The world would be a better place if people from other countries were more like Serbians”
t0 <- unname(table(mydf$better_world[mydf$cyrillic_second == 0]))
t1 <- unname(table(mydf$better_world[mydf$cyrillic_second == 1]))
mylabels <- c('Strongly Disagree', 'Disagree', "Doesn't Agree or Disagree", 'Agree', 'Strongly Agree')
par(mar = c(15, 5, 2, 2))
m1 <- barplot(prop.table(t1) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m1, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m0 <- barplot(prop.table(t0) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m10 <- barplot((prop.table(t1) - prop.table(t0)) * 100, ylim = c(-10, 10), ylab = 'Difference in % of Respondents',  border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -10, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)

######### "Generally speaking, Serbia is no better than most other countries"
t0 <- unname(table(mydf$most_places[mydf$cyrillic_second == 0]))
t1 <- unname(table(mydf$most_places[mydf$cyrillic_second == 1]))
mylabels <- c('Strongly Disagree', 'Disagree', "Doesn't Agree or Disagree", 'Agree', 'Strongly Agree')
par(mar = c(15, 5, 2, 2))
m1 <- barplot(prop.table(t1) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m1, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m0 <- barplot(prop.table(t0) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m10 <- barplot((prop.table(t1) - prop.table(t0)) * 100, ylim = c(-10, 10), ylab = 'Difference in % of Respondents',  border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -10, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)


######### "Generally, the more influence Serbia has on other nations, the better off they are"
t0 <- unname(table(mydf$influence[mydf$cyrillic_second == 0]))
t1 <- unname(table(mydf$influence[mydf$cyrillic_second == 1]))
mylabels <- c('Strongly Disagree', 'Disagree', "Doesn't Agree or Disagree", 'Agree', 'Strongly Agree')
par(mar = c(15, 5, 2, 2))
m1 <- barplot(prop.table(t1) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m1, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m0 <- barplot(prop.table(t0) * 100, ylim = c(0, 50), ylab = '% of Respondents', border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -3, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)
m10 <- barplot((prop.table(t1) - prop.table(t0)) * 100, ylim = c(-10, 10), ylab = 'Difference in % of Respondents',  border = FALSE, cex.axis = 1.75, cex.lab = 2)
text(mylabels, x = m0, 
     offset = -0.1, y = -10, srt = 45, xpd = TRUE, pos = 2, cex = 1.8)


############################################################### 
############ SI A.8 - H1: Models with Covariates ##############
###############################################################


################ Models with controls####################
#------------------ Unified DV ------------------------#
# Model for H1, with controls
m1a <- lm(Y ~ cyrillic_second  + cyrillic_start +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, with controls
m1b <- lm(better_world ~ cyrillic_second  + cyrillic_start +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, with controls
m1c <- lm(most_places ~ cyrillic_second  + cyrillic_start +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, with controls
m1d <- lm(influence ~ cyrillic_second  + cyrillic_start +
            gender + college + ageCat + married + unemployed + urban +
            other_languages + dif_lc_home + dif_lc_work + cyrillicImportance +
            serbian + sns + polKnowQ1 + polKnowQ2 + intPolitics, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     type = 'text',
                     se = list(se1a, se1b, se1c, se1d), 
                     covariate.labels = c('Treatment 1 (Survey Cyrillic)', 
                                          'Alphabet Ad (Cyrillic)', 
                                          "Gender = Man", 
                                          "Gender = Prefer not to share", 
                                          "Education = College", 
                                          "Age = 36-49", 
                                          "Age = 50-64", 
                                          "Age = 65+", 
                                          "Marital Status = Married", 
                                          "Unemployed", 
                                          "Urban", 
                                          "# of languages spoken = 1", 
                                          "# of languages spoken = 2", 
                                          "# of languages spoken > 2", 
                                          "Delta use of Cyrillic to Latin (at home)", 
                                          "Delta use of Cyrillic to Latin (at work)", 
                                          "Importance of Cyrillic",
                                          "Ethnicity = Serbian", 
                                          "Party Identification = SNS",
                                          "Knows the # of MPs", 
                                          "Knows how long is the presidential term", 
                                          "Interest in Politics"),
                     star.cutoffs = c(0.05, 0.01),
                     omit.table.layout = "n",	
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is bettern than most places", "Better World if Equal to Serbia"),
                     title = "Effect of Cyrillic Alphabet on Respondents' Own Attitutes, Including Covariates", 
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr'
)


############################################################################
################# SI A.9 - H1: Result by Ad's Alphabet #####################
############################################################################

# Create Dummies (start_treat1)
mydf$latin_latin <- as.numeric(mydf$latin_start == TRUE & mydf$latin_second == TRUE)
mydf$latin_cyrillic <- as.numeric(mydf$latin_start == TRUE & mydf$latin_second == FALSE)
mydf$cyrillic_latin <- as.numeric(mydf$latin_start == FALSE & mydf$latin_second == TRUE)
mydf$cyrillic_cyrillic <- as.numeric(mydf$latin_start == FALSE & mydf$latin_second == FALSE)

# Latin + Latin is the control

#------------------ Unified DV ------------------------#
# Model for H1, no controls
m1a <- lm(Y ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, no controls
m1b <- lm(better_world ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, no controls
m1c <- lm(most_places ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, no controls
m1d <- lm(influence ~ latin_cyrillic + cyrillic_latin + cyrillic_cyrillic, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))


# Table
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     covariate.labels = c('Ad = Latin, Alphabet = Cyrillic',
                                          'Ad = Cyrillic, Alphabet = Latin',
                                          'Ad = Cyrillic, Alphabet = Cyrillic'),
                     title = "Effect of Alphabet on Respondents' Own Attitudes",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01),
                     type = 'text'
)

#################################################################################################
############ SI A.10 - H1: Interaction Political Knowledge and Interest in Politics ##############
#################################################################################################

### Interest in Politics

#------------------ Unified DV ------------------------#
# Model for H1, interaction with interest in politics
m1a <- lm(Y ~ cyrillic_second * intPolitics, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, interaction with interest in politics
m1b <- lm(better_world ~ cyrillic_second * intPolitics, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, interaction with interest in politics
m1c <- lm(most_places ~ cyrillic_second * intPolitics, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, interaction with interest in politics
m1d <- lm(influence ~ cyrillic_second * intPolitics, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table 
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     type = 'text',
                     covariate.labels = c('Alphabet = Cyrillic',
                                          'Interest in Politics',
                                          'Alphabet = Cyrillic x Interest in Politics'),
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is better than most places", "Better world if equal to Serbia"),                      
                     title = "Effect of Alphabet on Respondents' Own Attitudes Conditional upon Interest in Politics",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01)
)

## Create plots
model_list <- list(m1a, m1b, m1c, m1d)
sapply(model_list, plot_interaction)

#### Political Knowledge

#------------------ Unified DV ------------------------#
# Model for H1, interaction with political knowledge
m1a <- lm(Y ~ cyrillic_second * polKnow, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, interaction with political knowledge
m1b <- lm(better_world ~ cyrillic_second * polKnow, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, interaction with political knowledge
m1c <- lm(most_places ~ cyrillic_second * polKnow, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, interaction with political knowledge
m1d <- lm(influence ~ cyrillic_second * polKnow, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table 
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     type = 'text',
                     covariate.labels = c('Alphabet = Cyrillic',
                                          'Political Knowledge',
                                          'Alphabet = Cyrillic x Political Knowledge'),
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is better than most places", "Better world if equal to Serbia"),                      
                     title = "Effect of Alphabet on Respondents' Own Attitudes Conditional upon Political Knowledge",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01)
)

## Create plots
model_list <- list(m1a, m1b, m1c, m1d)
sapply(model_list, plot_interaction, z_values = 0:2, z = 'cyrillic_second:polKnow', z_name = 'Political Knowledge')

#### Party ID

#------------------ Unified DV ------------------------#
# Model for H1, interaction with partisan preference
m1a <- lm(Y ~ cyrillic_second * sns, data = mydf)
se1a <- sqrt(diag(vcovHC(m1a, 'HC2')))

#------------------ better_world ------------------------#
# Model for H1, interaction with  partisan preference
m1b <- lm(better_world ~ cyrillic_second * sns, data = mydf)
se1b <- sqrt(diag(vcovHC(m1b, 'HC2')))

#------------------ most_places ------------------------#
# Model for H1, interaction with partisan preference
m1c <- lm(most_places ~ cyrillic_second * sns, data = mydf)
se1c <- sqrt(diag(vcovHC(m1c, 'HC2')))

#------------------ influence ------------------------#
# Model for H1, interaction with partisan preference
m1d <- lm(influence ~ cyrillic_second * sns, data = mydf)
se1d <- sqrt(diag(vcovHC(m1d, 'HC2')))

# Table 
stargazer::stargazer(m1a, m1b, m1c, m1d, 
                     se = list(se1a, se1b, se1c, se1d),
                     type = 'text',
                     covariate.labels = c('Alphabet = Cyrillic',
                                          'Party Identification = SNS',
                                          'Alphabet = Cyrillic x Party Identification = SNS'),
                     dep.var.labels = c('Nationalism index', "More Serbian is better", 
                                        "Serbia is better than most places", "Better world if equal to Serbia"),                     
                     title = "Effect of Alphabet on Respondents' Own Attitudes Conditional upon Party Identification = SNS",
                     single.row = TRUE,
                     no.space = TRUE ,
                     omit.stat = c("ser", 'f'),
                     style = 'apsr',
                     star.cutoffs = c(0.05, 0.01)
)

## Create plots
model_list <- list(m1a, m1b, m1c, m1d)
sapply(model_list, plot_interaction, z_values = 0:1, z = 'cyrillic_second:sns', z_name = 'Party Identification = SNS')
